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Abstract 

We present results from two different methods for the calculation of hadron spectra in 
QCD and SUSY QCD with large primary energies ^/s up to 10^^ GeV. The two methods 
considered are a Monte Carlo (MC) simulation and the evolution of fragmentation func- 
tions described by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations. 
We find that the pion, nucleon and all-hadron spectra calculated with the two methods 
agree well. The MC simulation is performed with new hadronization functions (in com- 
parison with our previous work), motivated by low-energy (y^ < Mz) data and DGLAP. 
The hadron spectra calculated with both sets of hadronization functions agree well, which 
indicates that our method for calculating the hadronization function works successfully. 
The small difference in the calculated hadron spectra characterizes the uncertainties of 
this method. We calculate also the spectra of photons, neutrinos and nucleons and com- 
pare them with other published results. The agreement is good for all x from ~ 10~^ 
up to X < 0.3. The consistency of the spectra calculated by different methods allows 
to consider the spectral shape as a signature of models with decays or annihilations of 
superheavy particles, such as topological defects or superheavy DM. The UHECR spectra 
from these sources are calculated. 
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1 Introduction 



Ultra-High Energy Cosmic Rays (UHECR) remain a puzzle in physics. First of all, eleven 
AGASA events with ii^ > 1 x 10^*^ eV contradict the Greisen-Zatsepin-Kuzmin (GZK) 
cutoff 0, although the HiRes data [3] are generally considered to be consistent with the 
GZK cutoff [1]. If the UHECR primaries are protons (see below) and if they propagate 
rectilinearly, as the claimed correlations with BL Lacs at lower energies (4— 8) x 10^^ eV jS] 
imply, then their sources must be seen in the direction of the highest energy events with 

~ (2 - 3) X lO^o eV detected by HiRes 0, Fly's Eye 6_ and AGASA 0. Indeed, 
the proton attenuation length at these energies is only 20 - 30 Mpc and the sources 
should have been seen in the direction of these particles, since such correlations exist at 
considerably lower energies. This implies that particles with E ~ lO^'^ eV may have a 
different origin as those with lower energies. 

Meanwhile, there is strong evidence that primary particles at lower energies, 1 X 10^^ < 

< (7 — 8) X 10^^ eV, are extragalactic protons, most probably from Active Galactic 
Nuclei (AGN). The different pieces of evidence include: (i) Extensive Air Shower (EAS) 
data confirm protons as primaries (ii) the dip seen with Xdof — ^■'^ ™ 

spectra of AGASA, HiRes, Fly's Eye and Yakutsk , is a signature of the propagation 
of UHE protons in the extragalactic space and (iii) the beginning of the GZK cutoff seen 
in the spectra of AGASA and HiRes ^2j . 

According to the correlations with BL Lacs found in Ref. [H] and the analysis of small- 
angle clustering jl3j . protons should propagate rectilinearly from the sources (AGN). 
However, if one excludes the correlations with BL Lacs from the analysis, the propagation 
of protons in very strong magnetic fields becomes also feasible |14l I15j . Nevertheless, 
also in this case, the lack of a nearby source in the direction of the highest energy events 
(e.g. at -E ~ 3 X lO^'' eV) remains a problem for reasonable field strengths 5 ~ 1 nC: 
the deflection angle, 9 ~ lait/rn = 3.7° B^q given by the attenuation length /att and the 
Larmor radius rn, is small and sources should be seen. 

Many ideas have been put forward aiming to explain the observed superGZK (E > 
(6 — 8) X 10^^ eV) events: strongly interacting neutrinos and new light hadrons jl7j 
as unabsorbed signal carriers, Z-bursts J^, Lorentz-invariance violation ^j. Topological 
Defects (TD) (see EHl for a review), and Superheavy Dark Matter (SHDM) (see EH for 
a review). 

The two last models listed above share a common feature: UHE particles are produced 
in the decay of superheavy (SH) particles or in their annihilation. In the case of TD they 
are unstable and in the case of SHDM long-lived particles. We shall call them collectively 
X particles. In the Z-burst model, the decay of a much lighter particle, the Z-boson, is 
involved, while the decay products are boosted by very large Lorentz factors. Annihilation 
takes place in the case of monopolonia 22 , necklaces |23) and SHDM particles within some 
special models |24j . As elementary particle physics is concerned, both processes proceed 
in a way similar to e~^e~ annihilation into hadrons: two or more off- mass-shell quarks and 
gluons are produced and they initiate QCD cascades. Finally the partons are hadronized 
at the confinement radius. Most of the hadrons in the final state are pions and thus the 
typical prediction of all these models is the dominance of photons at the highest energies 
£; > (6 - 8) X 10^^ eV. 

This prediction is questioned by the AGASA data at > 1 x 10^^ eV combined 
with the recent recalculation (221 of the muon number in photon-initiated EAS at the 
highest energies. Previously |271 128j . the muon content in photon- induced showers was 
found to be similar to that in proton-induced showers. According to the new calculations. 



the muon number is still large, but 5-10 times smaller than in hadronic showers. From 
eleven AGASA events at > 1 x 10^'^ eV, the muon content is measured in six and in 
two of them it is high. The muon content of the remaining four events is close to the 
one in photon-induced showers. MC simulations show that a primary flux comprised 
completely by photons is disfavored. We shall discuss further these data in the Conclusions. 

The spectrum of hadrons produced in the decay /annihilation of X particles is another 
signature of models with superheavy X particles. Several methods for the calculation of 
these spectra have been developed in the past several years. 

The mass of the decaying particle, Mx, or the energy of annihilation ^/s, is in the 
range 10^^ - 10^^ GeV. The existing QCD MC codes become numerically unstable at 
much smaller energies, e.g., at Mx ~ 10'^ GeV in the case of HERWIG. Moreover, the 
computing time increases rapidly going to larger energies. Nevertheless, one of the first 
spectrum calculation has been performed with the help of HERWIG for energies up to 
Mx ~ 10^^ GeV, and the computed spectra were extrapolated to Mx ~ 10^^ GeV PU] . 
Supersymmetry was not included in these calculations. 

Another option used in the first calculations is the limiting spectrum, an analytic 
method developed in Refs. [HOI- This method has been found to be very successful at 
LEP energies (see j31j and references therein). Two basic assumptions are involved in this 
method: (i) the beta function (3 describing the running of the QCD coupling is taken 
to be constant, i.e. as(/c^) oc l/ln(A;^/A^) for all transverse momenta k_\_, and (ii) the 
minimum virtuality Qq of partons, down to which the cascade develops perturbatively, 
is taken equal to the scale A. The high energy supersymmetric generalization of this 
solution has been obtained in Ref. j32j. Later, a comparison with the MC simulation I 33j 
showed that the limiting spectrum does not describe well hadron spectra in SUSY QCD 
and that the assumption (i) is mainly responsible for the discrepancies found. Indeed, 
changing the evolution of as{k\) an agreement of the spectra around the Gaussian peaks 
was obtained when in the SUSY QCD MC with (5 = const, was used and 
was frozen at small k\, which is a reasonable physical assumption, e.g., in the DGLAP 
method^. For more details see Ref. j33j . 

Monte Carlo simulations are the most physical approach for high energy calculations 
which allow to incorporate many important physical features as the presence of SUSY 
partons in the cascade and coherent branching [^Hl- The perturbative part of our MC 
simulation is similar to other existing MC programs and hence reliable. 

For the non-perturbative hadronization part a phenomenological approach is used in 
Ref. The fragmentation of parton i into a hadron h is expressed through perturbative 
fragmentation function of partons Dj(x, Mx) convoluted with hadronization functions 
fj'{x,Qo) at the scale Qo, 

D'l{x,Mx)= E / -Di{x/z,Mx,Qo)f^{z,Qo), (1) 

where the hadronization functions do not depend on the scale Mx (as notation for the 
scale in this paper we shall use also = Mx and t = ln(s/so))- This important prop- 
erty of hadronization functions allows us to calculate //'(x, Qo) from available LEP data, 
D^{x, Mx) at the scale Mx = Mz, and then to use it for the calculation of fragmentation 
functions D^{x,Mx) at any arbitrary scale Mx- Our approach reduces the computing 
time compared to usual MC simulations and allows the fast calculation of hadron spectra 

We are not able to find any place in our paper [321 where according to [33j we "acknowledge errors.' 
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for large Mx up to Mqut- The MC of Ref. has two versions: one for ordinary QCD 
and another one for SUSY QCD. 

The perturbative part of the MC simulation in Ref. [SSI includes standard features 
such as angular ordering, which provides the coherent branching and the correct Sudakov 
form factors, as well as SUSY partons. Taking into account SUSY partons results only 
in small corrections to the production of hadrons, and therefore a simplified spectrum of 
SUSY masses works with good accuracy. 

The weak influence of supersymmetry is explained by the decay of SUSY partons, 
when the scale of the perturbative cascade reaches the SUSY scale QsusY ~ ^ TeV^. 
Most of the energy of SUSY partons remains in the cascade in the form of energy of 
ordinary partons, left after the decay of SUSY partons. The qualitatively new effect 
caused by supersymmetry is the effective production of the Lightest Supersymmetric 
Particles (LSP), which could be neutralinos or gluinos. The spectra of neutralinos were 
calculated in jSS]- 

MC simulations allow also to calculate the characteristic feature of the QCD spectrum, 
the Gaussian peak, which is beyond the power of the next method we shall review. 

The fragmentation functions D^{x, Mx) at a high scale Mx can be calculated evolving 
them from a low scale, e.g. Mx = Mz, where they are known experimentally. This evo- 
lution is described by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation 

\'A7\ which can be written schematically as 

3 

where t = ln(s/so), <8> denotes the convolution f ®g = dx/xf{x)g{x/z), and Pij is 
the splitting function which describes the emission of parton j by parton i. Apart from 
the experimentally rather well determined quark fragmentation function Dg{x, Mz), also 
the gluon FF D^{x,Mz) is needed for the evolution of Eq. @. The gluon FF can be 
taken either from MC simulations or from fits to experimental data, in particular to the 
longitudinal polarized e'^e~ annihilation cross-section and three-jet events. 

The first application of this method for the calculation of hadron spectra from decaying 
superheavy particles has been made in Refs. |HH1 EH]) followed by Refs. 1^ El HH 
1431 144j . The most detailed calculations have been performed in Ref. [1^, where more 
than 30 different particles were allowed to be cascading and the mass spectrum of the 
SUSY particles was taken into account. Although at Mz, which is normally the initial 
scale in the DGLAP method, the fragmentation functions for supersymmetric partons are 
identically zero, they can be calculated at larger scales t: SUSY partons are produced 
above their mass threshold, when their splitting functions are included in Eq. In this 
work we prove that this method is correct. Also, the LSP spectrum can be computed 
within the DGLAP approach O . 

A problem of the DGLAP method are the fragmentation functions at small x = 2p/ ^/s, 
where p is the momentum. The DGLAP equations allow to evolve fragmentation functions 
known at x > Xmin from a starting scale ^/sQ to a higher scale ^/s. Therefore, Xmin of the 
FF at the starting scale ^/sq determines the x range accessible for all s. Requiring that 
perturbation theory can be used, pj_ > A ~ 0.25 GeV, leads for the starting scale Mz to 
X > 0.005. 

In Ref. |43| I44j the initial fragmentation functions are taken from Ref. [IH] at the scale 
Qo ~ 1.4 GeV and are extrapolated to very low x ~ 10~^, i.e. into the non-perturbative 
region. The formal DGLAP evolution between the scales Qq and Mz is described by 
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equations with as{s) not depending on x (see e.g. Eq. (jlOp in Section I^J. Surprisingly 
enough this method works weU: The fragmentation functions for x as low as 10~^ evolved 
in Ref. [431 144j to large Mx coincide with our MC simulation. 

Near the GUT scale all three gauge couplings (i = 1, 2, 3) have approximately 
equal numerical values. Naively one expects that at scales Mx close enough to the GUT 
scale all particles including e.g. leptons and electroweak (EW) gauge bosons are cascading 
like QCD partons. 

In fact, the influence of the masses of the EW gauge bosons on soft singularities has 
to be carefully studied. In the MC approach, the leading effect of the finite masses of the 
gauge boson can be implemented in a rather straightforward way [17j- Cascading of the 
longitudinal modes of the EW gauge bosons or of the Higgs bosons is a subleading effect. 
In Ref. |17], it was demonstrated that EW cascading occurs at Mx > 10^ GeV, even if only 
leptons and EW bosons are included in the consideration. The interactions with quarks 
and gluons mixes EW and QCD cascades. If, for example, a hypothetical X particle 
couples at tree level only to neutrinos, their further cascading results in the production 
of the QCD partons and thus of hadrons. Thus the production of neutrinos through 
decays of X particles with such couplings is constrained by the usual electromagnetic 
(e-m) cascade limit in the universe. EW cascading has been included in the calculations 
of Refs. |43| I44j . though in a formal way. 

In this paper we shall study the agreement of two methods: MC and DGLAP equations 
for the calculation of spectra produced in the decay or annihilation of superheavy particles. 
For this we shall calculate the spectra using the same assumptions in both methods. We 
shall also compare the results obtained by different groups and confront the calculated 
spectra with recent ones measured by UHECR experiments. 

The paper is organized as follows: in Section 2 we introduce the DGLAP equations 
and the technique used to solve them. In Section 3 we discuss the properties of the 
hadronization functions used in the MC simulation and obtain new low-energy motivated 
hadronization functions as extension of those used in [HSI . Then we compare the fragmen- 
tation functions calculated with the MC and the DGLAP equations in Section 4. Photon, 
neutrino and proton spectra, needed for UHECR calculations, are computed in Section 
5 and compared with the spectra obtained in Refs. |34 | I43 | IH| . Finally, in Section 6 we 
consider the consequences of our results for models of superheavy DM and Topological 
Defects. 

2 DGLAP equations in SUSY QCD 

If the X particle decays into partons i = u,u, . . . , g, the parton EEs D^{x,m\) can be 
defined as the probability of fragmentation of a parton i into a hadron h with momentum 
fraction x = 2p/AIx- The evolution of EEs with increasing scale s = (or t = Ins/so) 
is governed by the DGLAP equations, 

dtD^{x,t) = / -DHx/z,t)P,^^{z,t) . (3) 

Multiplying Eq. (j^j) by x and integrating it over x shows that the DGLAP equations 
conserve momentum, dt J Dj{x,t)xdx = 0, if J2j I Pi^j{x)xdx = 0. 

Since in the limit s ^ all quark flavors couple to gluons in the same way, the gluon 



FF mixes only with the flavor singlet FF of quarks, 

J flavors 



(4) 



where the summation goes over the number of active quark flavors nj involved in the 
process (nj increases with increasing scale t). 

The coupled evolution equation for the gluon FF and the quark singlet FF 
becomes then 



D'^{x,t) 



D^{x,t) 



Pqq{X)t) Pgq{x,t) 
2nfPqg{x,t) Pgg{x,t) 







Dhx,t) 



(5) 



where (8) denotes as usually the convolution [/ 17] (z) = dx/xf{x)g{x/ z) = 
JI dx/xf{x/z)g{x), and Pij = Pi^j. 

A formal extension of the standard DGLAP equations © to the SUSY case is straight- 
forward. Denoting squark and gluino by q and g, respectively, and considering the flavor 
singlet FF for squarks as it was discussed for quarks, we can write the SUSY DGLAP 
equations as 



Dhx,t) 
Dhx,t) 
V Dlix,t) J 



( Pqg{x,t) Pgq{x,t) Pgq{x,t) Pgq{x,t) \ ( D^{x,t) \ 

Dg{x,t) 



2nfPqg{x,t) Pgg{x,t) 2nfPqg{x,t) Pgg{x,t) 
Pqq{x,t) Pgq{x,t) Pqq{x,t) Pgq[x,t) 
_\ 2nfPqg{x,t) Pgg{x,t) 2nfPqg{x,t) Pgg{x,t) J 



D\{x,t) 
V Dlix.t) j 
(6) 

The exact form of the splitting functions Pij{x, t) is not known, but using the perturbative 
expansion of these quantities one has 



n=0 



27r 



2-jr 



(7) 



The LO splitting functions of SUSY QCD were derived in Ref. [48) and are given in their 
unregularized form Pij in Tabled 

In the case of the diagonal splitting functions, we have to distinguish between regu- 
larized splitting functions Pii{x) and unregularized ones Pii{x). These splitting functions 
have a probabilistic interpretation only for x < 1, since they contain a delta function 
contribution at a; = 1 accounting for losses. Moreover, they describe the emission of 
soft gluons for x — > 1 and contain therefore a pole of infrared type which needs to be 
regularized. Using momentum conservation. 



dz z [P(^\z) + + + = 1 



gq 



I' dz z [2n;PW(z) + PW(Z) + 2n;p(°)(z) + pf^{z)\ = ] 
£ dz z [P^fiz] + P^f(z) + P^fiz) + P^l\z)] = 1 

£ dz z [2nfP^fiz) + P^fiz) + 2nfP^f{z) + P^f {z)] = \ 
one obtains as formal expression for the regularized splitting functions 
Pf^{^) = P!i?{x) - 5(1 - x) j\z z [PW(Z) + PW(Z) + Pf^{z) + P(;)(z)] 



(8b) 



(8d) 



(9a) 



(x) = P W (x) - 5(1 - x) 1^ z [2n^pW (.) + Pjo) (z) + 2n;p(°) (z) + pj^^ (z)] (9b) 

(^) = 4°^ (^) - '^(l - ^) J/^ ^ [P^f (^) + Pgf (^) + + (^)] (9^) 

P§f (^) = 4°^ (^) - ^ (1 - ^) /rf^ - [2n/Pi°^ (^) + P^f (z) + 2n;P|) (z) + (z)] . (9d) 

Instead of calculating explicitly the expressions after the delta functions, we substitute 
Eqs. (|9aM9d|) directly into the DGLAP equations. In the case of ordinary QCD, the 
evolution of the singlet quark FF is then given by 




where denotes the usual step function. Since the two terms in the square bracket cancel 
each other for z = 1, the pole of Pgq\z) disappears. The same method can be used to 
replace the other diagonal splitting functions by their unregularized counter-parts. 
Energy conservation is automatically ensured by Eqs. H9aH9d|) . 

The DGLAP equations allow to evolve the FFs -Df (x) known at some scale sq to higher 
scales s. Since we are interested in a comparison of this method with our MC simulation, 
we use for the initial scale sq = -^f; the same as we have used in the MC to derive 
the hadronization function. Alternatively, we shall use for the initial scale of evolution 
also the very low value y/s^ 1 GeV as it will be explained in Section |31 

In the case of the initial scale Mz one has two initial FFs, D^{x,Mz) and D^{x,Mz). 
The former is known experimentally and the latter is calculated using our MC simula- 
tion. We shall describe now shortly some technical aspects connected with the numerical 
solution of the DGLAP equations. Starting from an initial set of FFs at given t and x, 
all FFs are evolved simultaneously with a 4th order Runge-Kutta algorithm with fixed 
step-size. At each Runge-Kutta step, the rhs of the DGLAP equations is evaluated with 
a Gaussian quadrature algorithm. Since this algorithm requires the knowledge of the FFs 
at X values different from the initially chosen ones, a polynomial interpolation algorithm 
is used to calculate the FFs. To avoid the 1/z singularities in the integrand, we evolve 
xD^{x,t) instead of Df(x,t). 

For the evolution of as{s) as a function of s we use the same method as in our MC 
simulation [HSj : we combine the thresholds of gluinos and all squarks in a single threshold 
at s = MgugY- The numerical value of MsusY is then fixed by requiring unification of 
coupling constants as in the minimal SUSY SU(5) model. This simplified treatment of the 
SUSY mass spectrum allows to compare the two methods using the same assumptions. 
Since moreover the results depend only weakly on MgusY; this simplification is physically 
reasonable. 

We shall finish this section with a remark on the connection between the more of- 
ten discussed space-like evolution of structure functions / and the time-like evolution of 
fragmentation functions D describing (SUSY) QCD cascades [IHI. Since Dj represents 
the fragmentation of the final parton j, while fi describes the distribution of the initial 
parton i, the matrix of splitting functions P has to be transposed going from one case to 
the other, as 

(Pi^i)-- rv ^ (Pi^i)--,- i-v • (11) 

'ij, space— like J ''j i, time— like ' 



On the other hand, a formal analytic continuation relates the splitting functions in 
both regions at leading order (LO): neglecting color factors, this relation is 



space— like 



time— like 



(12) 



Performing both transformations, the DGLAP equations are identical in the time- and 
space- like region at LO. 

Since in Refs. |34[ I45j only the transformation Ullj) has been performed, the splitting 
functions there (most notably for gluons and gluinos) are different from ours and from 
those in Refs. [431 144j . where both transformations and (|12p have been correctly 
used. 



3 Hadronization functions in the Monte Carlo 
simulation 

In this Section we discuss the general properties of the hadronization functions used in 
the MC simulation and their connection with the DGLAP method. 

In the MC simulation the hadronization functions are defined by Eq. . They can be 
determined from the FFs at lower scales, e.g at Mx = Mz, known from e~^e~ — > hadrons 
data. Namely, we take the measured all-hadron spectrum as FF Dg{x, Mz) as Ihs of 
Eq. and compute the hadronization functions at the rhs of this equation. 

However, hadronization functions cannot be calculated in an unique way using Eq. (jTj): 
Only if the FF function D^{x,Mz) were known precisely at an infinite number of points 
X, then the hadronization functions fj{x) could be calculated in principle in an unique 
way using e.g. the method of inversion. In practice the number of points x, where the 
FFs are measured is limited, and the experimental errors at some of them (most notably 
at X close to 1) are large. Another uncertainty arises if the flavor dependence of FFs is 
considered. The arbitrary choice of the minimal virtuality Qq of partons in the MC also 
introduces some additional error: Although the rhs of Eq. should not depend on Qq) 
in practice all MC simulations yield slightly different results for different Qq. 

Instead of using the inversion method, we assume a specific functional form of the 
hadronization function, characterized by a set of free parameters, and perform then a fit 
to the data. In this method additional uncertainties appear, because the functional from 
has to be specified a priori. As a result one can obtain different hadronization functions 
within the uncertainties discussed above. Deriving a new set of hadronization functions 
in this paper, we estimate as by-product the corresponding uncertainties in the calculated 
spectra. 

In Ref. we have used two hadronization functions — one for the quark flavor singlet 
and another for the gluons — motivated by the limiting spectrum. Each of them contains 
three free parameters and we have determined these parameters using the observed spec- 
trum of charged hadrons at y/s = Mz- Then we performed several tests to check our 
hadronization scheme (see Ref. |33j). Among them are the calculation of hadron spectra 
at two other scales, ^/s = 58 GeV and ^/s = 133 GeV, where measurements are available. 
We also mention here one particular test based on the limiting spectrum. We run our 
ordinary QCD MC with the hadronization functions fixed as above using the same set 
of assumptions which are used in the derivation of the limiting spectrum: we fixed the 
number of flavors to nj = 3 (or to nj = 6) and the running of Ug was taken exactly as in 



the limiting spectrum method. We obtained an exceUent agreement at ah scales Mx and 
for all X except those close to one, where the limiting spectrum is invalid. 

The FFs in these calculations, D^{x,Mx) and Dj{x,Mx), conserve momentum ex- 
actly and with good accuracy, respectively. 

In this Section we shall use another choice of hadronization functions, imposing low- 
energy data and some additional physical interpretation motivated by the DGLAP equa- 
tions, and demonstrate that the calculated spectra for different Mx agree well with our 
previous calculations. 

We shall begin with the general properties of fragmentation and hadronization func- 
tions, which we impose on those used in our calculations. 
(i) Momentum conservation of hadrons, 

J2h lo ^ii^' s)xdx = 1, and partons J2j lo ^1 {^^ s)xdx = 1, results in the normalization 
of the hadronization function as ^hlo fi{x)xdx = 1. 
The proof follows from Eq. (^. Indeed, 



1 = ^ /' s)xdx = Y.T.t /' 

^ Jo ^ . JQ Jx Z Z 



s)fHz). (13) 



Changing the order of integration, using the variable x' = x/z and introducing the 
new function ^ 

C/(s) = / D{{x',s)x'dx', (14) 

one obtains 







1 = EE [' dzzf^{z)C{{s), (15) 
h j 

which must hold for arbitrary s. Consistency with the condition J2j Cj (s) = 1 leads to 

[\zzft{z) = l (16) 

as the only solution. 

Similarly, one can prove the reversed statement: 
(ii) The momentum conservation of partons in the perturbatively calculated cascade 

Jq Dj{x, s)xdx = 1, with the hadronization function normalized as in Eg. Illb\) results 
in momentum conservation of hadrons. 

Now we add an assumption about hadronization functions which links MC and DGLAP 
methods, namely we assume that the hadronization function is a FF function D^{x, ^/s) 
extrapolated to the very low scale ^/s = Qq. Then the FF L>f(x,^) can be computed 
with the help of DGLAP evolution equations using the hadronization function //^(x, Qq) as 
the initial FF. Under this assumption we should choose the hadronization function as the 
fragmentation function obtained for low energy e~^e~ annihilation and/or ep scattering. 
Namely, we take it following Ref. [^0] as 

f^{x, Qo) = Nix'^^ {x + xi)-'^ (1 - xf^ , (17) 

where the index i runs through quark singlet q and gluon g, and the index h refers to 
the total hadron spectra. We fix the value of Qq in our MC simulation as in Ref. [SHI 
to Qo = 0.625 GeV^. As in the parameters are found performing a fit of the 
LEP spectrum at ^/s = Mz by D^{x,Mz) calculated from Eq. The hadronization 
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functions are shown in Fig. ^ Their main features agree with those of fg has its 
maximum at x ~ 0.1, while fq peaks close to x = 1. As a difference one can see that the 
new fit function is chosen to ensure fq{x = 1) = 0. 

In Figs, mini and 01 the spectra of hadrons D^{x,y/s) calculated for the scales y/s = 
58 GeV, 91.2 GeV, and 133 GeV are compared with observations. (The measurements 
are for charged hadrons and rescaled by us to the total hadron spectra.) In Figs. Eland 
ini these spectra computed at the scales y/s = 1 x lO^'^ GeV and 1 x 10^^ GeV using the 
old [HHl and new hadronization functions are shown. The agreement is good, as it has 
been expected. The ratio between the new and old spectra, shown in Fig. [3 illustrates 
the uncertainties in our MC simulations connected with different choices of hadronization 
functions. Figures El and [7| show that these uncertainties affect only the high energy 
part of the spectrum at x > 0.1. 

We described above the all-hadron spectra D^{x, Mx)- Similar calculations have been 
performed by us for charged pions and nucleons, using the experimental data at ^/s = Mz 
from the ALEPH and OPAL cohaborations |53j . 

As a test of the interpretation of the MC hadronization function /j^(x,(5o) as low 
energy limit of the FF L>f (x, ^/s), we have evolved the hadronization function given by 
Eq. ((T7|) from the scale Qq = 0.625 GeV^ to the higher scales using our QCD DGLAP code 
described in Section j^l In principle, this procedure should not work well for very small 
X = 2p/Qo, where one is physically in the non-perturbative regime (see Introduction). 
But using the DGLAP equations which prescribe s as argument of ag, we can perform 
formally the evolution calculations, as it was done in Refs. |46| I5()j and j43| I44j . 

In Fig. |Hland|nithe evolved hadronization functions at the scale Mz are compared with 
FFs L>f(x,Mz) calculated with the MC. We find good agreement for Dg but much worse 
agreement for D^ . We think that a reason for this failure is the use of the quark flavor 
singlet FF at low scales. Indeed, in the successful evolution of the Refs. [SlISllMlEni, the 
FFs for different quarks have different parameterizations. For scales above Mz, the flavor 
singlet FF becomes a good approximation, and, indeed, the evolution of hadronization 
functions from the scale Qq to scale Mgut results in very good agreement with MC 
simulation (see Fig. I10|) . In this case, the agreement is equally good for both and D^. 

4 Comparison of DGLAP and MC hadron spec- 
tra 

In this Section we shall compare the hadron spectra computed by the two methods dis- 
cussed above, MC and DGLAP. 

We shall begin with a remark concerning the smallest and largest x values of practical 
interest. Coherent branching produces the so-called Gaussian peak in multiplicity with 
maximum at x^ax ~ (Qo/Mx)^'^ ~ 2 x 10~^° for Mx ~ Mqut- This approximate 
estimate coincides well with the value found in the MC simulation [^. Note that the 
DGLAP method is not valid for those x values where coherence plays an essential role, 
namely at x < 4 x 10"^ (see Fig. 8 in Ref. 23])- The lowest value Xmin relevant for 
physical applications is much higher than Xmax- ^min ~ 2Eohs/Mx ~ 2 X 10-6 for the 
same Mx ~ Mqut and for -Eobs > 1 x lO^'^ GeV. On the other side, the maximal observed 
energy ~ 3 x lO^'' eV and the minimal Mx of interest, Mx > 10^^ GeV, restrict x to 
values much smaller than 0.6. Thus the value x = 0.6 can be considered as the largest x 
of interest for existing experimental data. 
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In Fig. ^Jwe present a comparison of D^{x, Mx) with i = q and Mx = 1 x 10^^ GeV 
calculated for ordinary QCD with the MC and DGLAP method. The initial scale in the 
DGLAP method is taken as ^/s = Mz, but the initial scale ^/s = Qq gives practically an 
identical spectrum. The spectra for i = g agree equally well, as well as the spectra for 
other high scales Mx- 

One can see that the MC and DGLAP spectra slightly differ at very low x and have 
a more pronounced disagreement at large values of x. The discrepancy at low x is due 
to coherent branching, and it starts at the value of x estimated above (x < 4 x 10~^). 
At large x, the calculations by both methods suffer from uncertainties, particularly the 
MC simulation. In this region the results are sensitive to the details of the hadronization 
scheme (see, e.g., the problem of HERWIG [SSI with the overproduction of protons at 
large x and the dependence on the choice of the hadronization functions in our MC as 
illustrated by Figs. 0and|ni). One can add to this problem large uncertainties in the 
measured FFs at Mz and x ~ 1 and also the theoretical uncertainties connected with 
the models of X particle decay (unknown number and types of the initial partons^ and 
unknown matrix element of the X particle decay). 

However, as a whole, Fig. II II demonstrates good agreement between MC and DGLAP 
methods. In Fig. ^1 we present the ratio of FFs calculated with MC and DGLAP for 
ordinary QCD. At x < 1 x 10"^ the MC FF is noticeably smaller than the DGLAP 
FF because coherence effects suppress branchings, an effect not included in the DGLAP 
method. At large x > 0.1, the discrepancy becomes greater than 20% due to the reasons 
explained above. For 2.5 x 10~^ ^ x < 0.1, the agreement between these two FFs is better 
than 20%. 

Let us now come over to SUSY QCD. 

In Fig. iniwe plot the FFs D^{x,Mx) calculated by MC and DGLAP methods for 
Mx = 1 X 10^*^ GeV and i = q. In the DGLAP method the SUSY FFs have been evolved 
from the ones obtained with the SUSY MC at the scale ^/s = lOMsusY ^ 10 TeV. One can 
see the good agreement between DGLAP (solid curve) and MC (dotted curve). This good 
agreement holds also for other (lower) scales Mx and for other initial partons i = g,g,q. 
The ratio of FFs for SUSY MC and SUSY DGLAP is shown in Fig. HHJ 

When one does not have the initial SUSY FFs from a MC simulation, the question 
arises how to proceed. As was first suggested in Ref. IHEll^) the initial FFs can be 
taken as the ones for ordinary QCD at the low scale ^/s = Mz, while the production 
of SUSY partons is included in the splitting functions assuming threshold behavior at 
-^SUSY ~ 1 TeV. We can check this assumption computing the SUSY FF in both ways. 
In Fig. HSlwe present the SUSY FFs D^{x,Mx) ior i = q and Mx = 1 x 10^^ GeV, 
evolved from the initial scale ^/s = Mz (dashed curve). The good agreement between the 
two DGLAP curves proves the validity of the assumption made above. 



5 Photon, neutrino and nucleon spectra 

The spectra of photons, neutrinos and nucleons produced by the decay of superheavy 
particles are of practical interest in high energy astrophysics. These spectra Df{x,Mx) 

particles decay in many models due to nonperturbative interactions, and the number of primary partons 
is therefore model dependent. Since at large x the number of cascade generations is small, the FFs are for x — > 1 
rather sensitive to the initial multiplicity. In contrast, at small x when the number of cascade generations is 
large, differences in the initial multiplicity become inessential. 
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with a = 'y,i',N can be also considered as FFs. Because the dependence on the type i of 
the primary parton is weak, we shah omit the index i from now on, keeping a as subscript. 

TiU now we concentrated our discussion on the total number of hadrons (a = h) 
described by the FF Dh{x, Mx), but in fact we have performed similar calculations sep- 
arately for charged pions and protons+antiprotons. The procedure of the calculations 
is identical to that described in Section [21 for the DGLAP method and in Section |31 for 
the MC. For charged pions and protons+antiprotons we used experimental data from 
Refs. [SHI- Below we shall present results of our SUSY MC simulations in terms of FFs 
for all pions Dt^, all nucleons D]\f and all hadrons Dh- We introduce the ratios £n{x) and 
as 

Dn(,x) = eN{x)Dh{x) , D^{x) = e^{x)Dh{x) . (18) 

The spectra of pions and nucleons at large Mx have approximately the same shape as 
the hadron spectra, and one can use in this case e-,^ = 0.73 it 0.03 and en = 0.12 it 0.02, 
taking into account the errors in the experimental data [511 1521 15^ . In Fig. E]the ratios 
eAr(x) and eT^{x) are plotted as functions of x for different values of Mx- Note the peculiar 
dependence of e^ix) for small Mx- The smallness of e^ix) at small x and small Mx is 
caused formally by the smallness of the observed vr/p ratio at Mx = Mz, where we fit our 
FFs. Physically it is due to the combined effect of coherence and the mass difference of 
nucleons and pions. Indeed, the x values at small Mx where £Nix) is particularly small 
belong to the region below the Gaussian peak {xm ~ (A/Mx)^'^), where coherence effects 
suppress branchings particularly strong. 

We can calculate now the spectra of photons and neutrinos produced by the decays 
of pions neglecting the small contribution (0.15 it 0.04) of K, D, A and other particles. 
Including these particles affects stronger neutrinos than photons, which are the main topic 
of this Section. For the pion spectrum we shall use D-,^{x) = s-,r{x)Dh{x). 

The normalized photon spectrum from the decay of one X particle at rest is given by 

D,{x) = '^ f'^eAy)DUy). (19) 
3 Jx y 

The total neutrino spectrum from decays of charged pions and muons can be repre- 
sented as [Hi] 

D,{x) = d:;^'^'^-{x) + D(;----=^(x) + z)^;f '^-'^'^^(x) , (20) 

where for pion decay 

d:;;>"'^{x) = \Rf ^ eAy)D,{y) (21) 
3 JxR y 

with 

R={l-ml/ml)-\ (22) 

and for muon decay 

2 (ill r^/fy 

^M-^M^ee(^) = -e.R / ^g.(y) / - e.(^) Dh{z) (23) 
<j -I X y -^^/y ^ 

with 

5 4 

9i = - - + -y^ for v^,^^ and = 2(1 - 3y^ + 2y^) for z/e,i^e- (24) 

The spectra are presented in Fig. for different masses Mx- We shall compare 
our photon spectra with those calculated by the DGLAP method in Refs. |431 144j and 
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|34j . The photon spectrum is most interesting to compare, because it is straightforwardly 
related to the hadron spectrum which is the basic physical quantity. Moreover, the photon 
spectrum is the dominant component of radiation produced by superheavy particles. 

To be precise, we compare the FF D'^{x,Mx) at Mx = 1 x 10^^ GeV. Figure [TBI 
demonstrates good agreement between our spectrum and those from Refs. [ISl IB] at 
X < 0.3. As was it mentioned above, the disagreement at large x is not surprising. Apart 
from Dg{x, Mz) taken directly from the experiments, both calculations use the much 
more uncertain D^{x,Q'^). In our case, D^{x,Q'^) is taken from our MC simulation |33j . 
in the case of Ref. |43[ I44j from the fit performed in Ref. j4tjj . In both cases, rather 
large uncertainties exist at large x (e.g. see Fig. 5 from Ref. [IHl)- In both calculations 
the decay of the X particle into two partons is considered, but in fact in many models 
only many-parton decays exist. This and the unknown type of the initial partons add a 
common theoretical uncertainty to both calculations (see Section 0] for discussion) 

The photon spectrum of Ref. [HI] shows some deviation at x < 0.3 from both spectra 
discussed above. To find the reason we performed the same calculations using splitting 
functions according the prescription of [31] (see Section I^J and obtained indeed some 
excess at small x and good agreement at large x, as one observes in Fig. 1161 Nevertheless, 
the agreement between the three curves as presented in Fig. ^Jis good. 

In Fig. El we present also the comparison of our proton spectrum, computed with the 
SUSY QCD MC, with that of Refs. [13 EH and IHl , which shows good agreement. 

6 UHECR from Superheavy DM and Topological 
Defects 

As follows from Sectional the accuracy of spectrum calculations has reached such a level 
that one can consider the spectral shape cLS cl SI gnature of the model. The predicted 
spectrum is approximately tx dE/E^'^ in the region of x at interest, and it is considerably 
steeper than the QCD MLLA spectrum used in the end of 90s. The generation spectra 
for nucleons, neutrinos and gammas are shown Fig. 1151 

Another interesting feature of these new calculations is a decrease of the ratio of 
photons to nucleons, in the generation spectrum. This ratio is presented in Fig. 

1181 for Mx = 1 X 10^^ GeV by a solid curve together with a band of uncertainties given 
by the two dashed curves. At ~ 1 x 10~^ this ratio is characterized by a value of 2 
- 3 only. The decrease of the "f/N ratio is caused by a decrease of the number of pions 
in the new calculations and by an increase of the number of nucleons. This result has 
an important impact for SHDM and topological defect models because the fraction of 
nucleons in the primary radiation increases. However, in both models photons dominate 
(i.e. their fraction becomes > 50%) at > (7 — 8) x 10^^ eV (see below). 

In this Section we shall consider two applications: superheavy dark matter (SHDM) 
and topological defects (TD). 

UHECRs from SHDMs have been suggested in Refs. |55| l5Hj and further studied in 
Ref. I^ni . Production of SHDM particles naturally occurs in the time- varying gravitational 
field of the expanding universe at the post-inflationary stage |57j . 

The relic density of these particles is mainly determined (at fixed reheating tempera- 
ture and inflaton mass) by their mass Mx- The range of practical interest is (3 — 10) x 10^'^ 
GeV, at larger masses the SHDM is a subdominant component of the DM. 
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SHDM is accumulated in the Galactic halo with the overdensity 
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(25) 



where p^^^ « 0.3 GeV/cm^, p„ = 1.88 x IQ-'^^h'^ g/cm^ and 1^cdm/i^ = 0.135 P]. With 
these numbers, 6 ~ 2.1 x 10^. Because of this large overdensity, UHECRs from SHDM 
have no GZK cutoff 

Clumpiness of SHDM in the halo can provide the observed small-angle clustering 1^ . 

The quotient rx = (to /tx) of relic abundance Qx and lifetime tx of the X particle 
is fixed by the observed UHECR flux as rx ~ 10~^^. The numerical value of rx is 
theoretically calculable as soon as a specific particle physics and cosmological model is 
fixed. In the most interesting case of gravitational production of X particles, their present 
abundance is determined by their mass Mx and the reheating temperature Tr. The 
life-time of the X particles is on the other hand fixed by choosing a specific particle 
physics model. As shown in Ref. there exist many models in which SH particles 
can be quasi-stable with lifetime tx 3> 10^" yr. The measurement of the UHECR flux, 
and thereby of rx, selects from the three-dimensional parameter space {Mx,Tji,tx) a 
two-dimensional subspace compatible with the SHDM hypothesis. Such a determination 
of a priori free parameters from experimental data has nothing to do with fine-tuning^, 
a reproach sometimes used against SHDM: Choosing a single value of tx from the wide 
range of theoretically allowed values just reflects the present state of theoretical ignorance. 

In Fig. I19| the spectra of UHE photons, neutrinos and protons from the decays of 
SHDM particles with Mx = 1 x 10^"^ GeV in the Galactic halo are presented. We have 
performed also a fit to the AGASA data using the photon flux from the SHDM model and 
the proton flux from uniformly disributed astrophysical sources. For the latter we have 
used the non-evolutionary model of Ref. jJOl- The photon flux is normalized to provide 
the best fit to the AGASA data at ^ > 4 x 10^*^ eV. The fits are shown in Fig. EOl with 
X^/d.o.f. indicated there. 

One can see from the fits in Fig. \2(J\. that the SHDM model with the new spectra can 
explain only the excess of AGASA events at E > 1 x W^^ eV: depending on the SHDM 
spectrum normalization and the details of the calculations for the extragalactic protons, 
the flux from SHDM decays becomes dominant only above (6 — 8) x 10^^ eV. 

Topological Defects (for a review see 63 ) can naturally produce UHE particles. The 
pioneering observation of this possibility has been made in Ref. |64j . 

The following TD have been discussed as potential sources of UHE particles: su- 
perconducting strings, ordinary strings, monopolonium (bound monopole-antimonopole 
pair), monopolonia (monopole-antimonopole pairs connected by a string), networks of 
monopoles connected by strings, vortons and necklaces (see Ref. [20] for a review and 
references). 

Monopolonia and vortons are clustering in the Galactic halo and their observational 
signatures for UHECR are identical to SHDM. However, as has been demonstrated in 
Ref. inS] ) the friction of monopolonia in cosmic plasma results in monopolonium lifetime 
much shorter than the age of the universe. 

Of all other TD which are not clustering in the Galactic halo, the most favorable 
for UHECR are necklaces. Their main phenomenological advantage is a small separation 
which ensures the arrival of highest energy particles to our Galaxy. 

•^In elementary-particle physics, fine-tuning is understood as tuning a quantity B to another quantity A with 
such precision that the predicted value m ~ A — B is many orders of magnitude smaller than A and B. 
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We shall calculate here the flux of UHECR from necklaces following the works [231 162j . 

Necklaces are hybrid TD produced in the symmetry breaking pattern G — > HxU{l) 
H X Z2. At the first symmetry breaking monopoles are produced, at the second one each 
(anti-) monopole get attached to two strings. This system resembles ordinary cosmic 
strings with monopoles playing the role of beads. Necklaces exist as the long strings and 
loops. 

The symmetry breaking scales of the two phase transitions, r/^ and r]s, are the main 
parameters of the necklaces. They determine the monopole mass, m ~ ATrrjm/e^ and the 
mass of the string per unit length /i ~ 27r7/g. The evolution of necklaces is governed by 
the ratio r ~ m/^d, where d is the average separation of a monopole and antimonopole 
along the string. As it is argued in Ref. j^S]) necklaces evolve towards configuration with 
r ^ 1. Monopoles and antimonopoles trapped in the necklaces inevitably annihilate in 
the end, producing heavy Higgs and gauge bosons {X particles) and then hadrons. The 
rate of X particles production in the universe can be estimated as |23j 

where t is the cosmological time. 

The photons and electrons from pion decays initiate e-m cascades and the cascade 
energy density can be calculated as 

1 2 dt 1 _ 3 2 A* 



where z is the redshift and ~ 1 is the fraction of the total energy release transferred 
to the cascade. 

The parameters of the necklace model for UHECR are restricted by the EGRET 
observations [HHI of the diffuse gamma-ray flux. This flux is produced by UHE energy 
electrons and photons from necklaces due to e-m cascades initiated in collisions with CMB 
photons. In the range of the EGRET observations, 10^ — 10^ MeV, the predicted spectrum 
is oc E~°' with a = 2 ^BT. The EGRET observations determined the spectral index as 
a = 2.10 lb 0.03 and the energy density of radiation as Wobs ~ 4 x 10~^ eV/cm^. The 
cascade limit consists in the bound Wcas < "^obs- 

According to the recent calculations of Ref. [HH] , the Galactic contribution of gamma 
rays to the EGRET observations is larger than estimated earlier, and the extragalactic 
gamma-ray spectrum is not described by a power-law with a = 2.1. In this case, the limit 
on the cascade radiation with a = 2 is more restrictive and is given by 

Wcas < 2 X 10"^eV/cm^ (28) 

We shall use this limit in further estimates. Using Eq. (|27jl with f^ = l and to = 13.7 Gyr 
ISHl we obtain from Eq. ^ r^fi < 8.9 x 10^"^ GeV^. 

The important and unique feature of this TD is the small separation D between 
necklaces. It is given Since r^^u is limited by e-m cascade radiation, 

Eq. (|27() . we can obtain a lower limit on the separation between necklaces as 



^'^ to > 10(^/10^ GeV2) V4 . (29) 



^ , .1/4 



This small distance is an unique property of necklaces allowing the unabsorbed arrival 
of particles with the highest energies. 
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The fluxes of UHECR from necklaces are shown in Fig |21l (the details of these cal- 
culations will be published in a forthcoming paper with P. Blasi [69 J. We used in the 
calculations r^// = 4.7 x 10^^ GeV^ which corresponds to uicas — 

1.1 X 10"^ eV/cm^, 

i.e. twice less than allowed by the bound The mass of the X particles produced by 
monopole-antimonopole annihilations is taken as Mx = 1 x 10^'* GeV. 

From Fig. [^one can see that in contrast to previous calculations j^, the necklace 
model for UHECR can explain only the highest energy part of the spectrum, with the 
AGASA excess somewhat above the prediction. This is the direct consequence of the 
new spectrum of particles in X decays obtained in this work. Thus UHE particles from 
necklaces can serve only as an additional component in the observed UHECR flux. 

7 Conclusions 

In this paper we have compared the MC and DGLAP methods for the calculation of 
hadron spectra produced by the decay (or annihilation) of superheavy X particles with 
masses up to Mgut ~ 1 x 10^^ GeV. We found an excellent agreement of these two 
methods. 

We have further elaborated the MC simulation of Ref. including the low-energy 
motivated hadronization functions //^(x, Qo) with the properties of fragmentation function 
D^{x,Qo) at low scale Qo ~ 1 GeV. Though the new hadronization functions are some- 
what different from the old ones, all new spectra computed for different initial partons i 
and different scales Mx agree very well with the old spectra, as it should be (see Section 
OJ. The small differences in the spectra illustrate the uncertainties involved in the ex- 
traction of hadronization functions from experimental data. However, these uncertainties 
affect the spectra only at x > 0.2. 

The calculations have been performed both for ordinary QCD and SUSY QCD. The 
inclusion of SUSY partons in the development of the cascade results only in small correc- 
tions, and it justifies our computation scheme with a single mass scale MsusY- 

In comparison to the DGLAP method, the MC simulation has the advantage of includ- 
ing coherent branching. It allows reliable calculations at very small x. The Gaussian peak, 
the signature of the QCD spectrum, cannot be obtained using the DGLAP equations. 

We have calculated the all-hadron spectra, as well as spectra of charged pions and 
nucleons, using the DGLAP equations. Two different methods (both based on our MC 
simulation) have been used. 

In the first one we have used the initial fragmentation functions, D^{x, Mz) and 
Dg{x, Mz), the former determined from hadron spectrum measured from e^e^-annihilation 
at ^/s = Mz and the latter calculated by MC. Then we evolved these fragmentation func- 
tions to higher scales Mx with the help of the DGLAP equations. In the second method 
we have used the hadronization functions as initial fragmentation functions, and evolved 
them from the scale Qo ~ 1 GeV to Mx- In the first method the spectra calculated 
at different scales and for different initial partons are in excellent agreement between 
themselves and with MC at x < 0.3 (see Figs. [TTl [T^ . 

In the second method the agreement is equally good at scales y/s ^ Mz- This also 
shows that hadronization functions can be seen as fragmentation functions at low scale 
Qo ~ 1 GeV. At small scales the quark singlet FF is calculated less precisely. 

The disagreement at x ~ 1 is explained by the uncertainties in the calculations using 
DGLAP and MC at large x, especially in the latter. Using the hadronic fragmentation 



1 c 



functions we have calculated the spectra of photons (produced by decays of neutral pions), 
neutrinos (from charged pions) and nucleons. 

Our nucleon spectrum agrees well with that of Refs. [IS] smd |44j . 

We compared also our spectrum of photons with the calculations of Ref. 
The comparison of the photon spectra is interesting, because of physical reasons (photons 
can be observable particles), and because the photon spectra are connected directly with 
the hadron spectra. 

The spectra are in good agreement (see Section El for the detailed discussion). The 
disagreement at the largest x is not of great practical interest because of the model 
dependent prediction of the spectrum. Indeed, because of the non-perturbative character 
of the decay, many-parton decays of X particles can dominate over the two parton decay 
considered (see Section [IJ) . Moreover, the range x > 0.3 corresponds to too high energies 
^ > 3 X 10^2 at the masses of superheavy particles at interest Mx > 1 x 10^^ GeV. 

We conclude that all calculations are in a good agreement especially at small x and the 
predicted shape of the generation spectrum (oc dE/E^'^) can be considered as a signature 
of models with decaying (annihilating) superheavy particles. 

The predicted spectrum of SHDM model cannot fit the observed UHECR spectrum at 
1 X lO^s eV < £" < (6 - 8) X 10^^ eV (see FigUHl). Only events at ^ > (6 - 8) x lO^^ gy, 
and most notably the AGASA excess at these energies, can be explained in this model. 
The robust prediction of this model is photon dominance. In present calculations this 
excess diminishes to ^/N ~ 2 — 3 (see Fig. I18|l . 

According to the recent calculation of Ref. , the muon content of photon induced 
EAS at -E > 1 X 10^*^ eV is high, but lower by a factor 5-10 than in hadronic showers. 
The muon content of EAS at > 1 x lO^'' eV has been recently measured in AGASA j25j . 
The measured value is the muon density at the distance 1000 m from the shower core, 
p^(lOOO). From 11 events at > 1 x 10^*^ eV the muon density was measured in 6. In 
two of them with energies about 1 x lO^'' eV, is almost twice higher than predicted for 
gamma-induced EAS. Taking into account the contribution of extragalactic protons at this 
energy (see Ref. |70| for an analysis), the ratio 7/p predicted by the SHDM model is 1.2 - 
1.4. It is lower than the upper limit 7/p < 2 obtained by AGASA at = 3 x 10^^ eV on 
the basis of a much larger statistics. The muon content of the remaining 4 EAS marginally 
agrees with that predicted for gamma-induced showers. The contribution of extragalactic 
protons for these events is negligible, and the fraction of protons in the total flux can be 
estimated as 0.25 < p/tot < 0.33. This fraction gives a considerable contribution to the 
probability of observing 4 showers with slightly increased muon content. Not excluding 
the SHDM model, the AGASA events give no evidence in favor of it. 

The simultanous observation of UHECR events in fluorescent light and with water 
Cherenkov detectors has a great potential to distinguish between photon and proton 
induced EAS. An anisotropy towards the direction of the Galactic Center is another 
signature of the SHDM model. Both kinds of informations from Auger |71j will be crucial 
for the SHDM model and other top-down scenarios. 

Topological defect models are another case when short-lived superheavy particle decays 
can produce UHECR. In Fig. |^ the spectra from necklaces are presented. One can see 
that at -E > 1 X 10^*^ eV photons dominate, and the discussion in the previous paragraph 
applies here too. In contrast to previous calculations [22]) the agreement with observations 
is worse: necklaces can explain only the highest energy part of the spectrum in Fig. I2H 
with the AGASA excess somewhat above the prediction. In the other energy ranges, UHE 
particles from necklaces can provide only a subdominant component. Other TDs suffer 
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even more problems (see Ref. |62|). 
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Table 1: Splitting functions for QCD and SUSY QCD from |^ {Cp = 4/3 and Nc = 3) 
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gluon 
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pW(x) = 2iVc + i + x(l - x) - 2 
P^^\x) = I [2x{l - x)] 
P^l\x) = 2Nc[x' + {l-xy] 


squark 


gluino 


P|)(x)=C^[2^] 
P^f{x) = Cp 


Pmi^) = i^ 

P^fi-) = Nc[^] 
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Figure 1: Hadronization functions for quarks fq{x, Qq) and gluons fg{x, Qq) obtained by fitting 
experimental data at ^/s = 91.2 GeV and using Ql = 0.625 GeV^. 




/ = ln(l/x) 

Figure 2: Comparison of the experimental data at y/s = 91.2 GeV and FF D^{x, Mz) computed 
by MC with hadronization functions as shown in Fig^ 
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Figure 3: Comparison of MC computed FF D^[x, ■\/s) with experimental data at ^/s — 133 
GeV. 
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Figure 5: Comparison of spectra calculated with the old |33] and new hadronization functions: 
FFs from ordinary QCD MC with quark as a primary parton are displayed for the scale ^/s = 
1 X 10^° GeV for the new (solid line) and old (dashed line) hadronization functions. 
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Figure 6: The same as in Fig. |31for the scale y/s = 1 x 10^^ GeV. 
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Figure 7: Ratio of the spectra calculated with the old and new hadronization functions for the 
scale ^ = 1 X 10^° GeV (solid line) and i/s = 1 x 10^^ GeV (dashed hne). For discussion of 
uncertainties at large x > 0.1 see the text. 




Figure 8: DGLAP FF with hadronization ftinctions as initial fragmentation functions at the 
scale Qq. The FF Dg{x, ^/s) is DGLAP evolved from the scale ^/s = Qq to the scale Mz 
(dashed curve) and plotted in comparison with MC calculated Dg(x, Mz) (full curve). 
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Figure 10: SUSY DGLAP FF with hadronization functions as initial fragmentation functions 
at the scale Qq. The FF -D^ (x, Mx) is SUSY DGLAP evolved from the scale y/s = Qo to scale 
Mx = 10^^ GeV (dashed curve) and plotted in comparison with MC calculated D'^{x, Mz) (full 
curve). 
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Figure 11: Comparison of the DGLAP FF (solid line) and MC FF (dashed hne) for ordinary 
QCD with Mx — 1 x 10^^ GeV and with quark as a primary parton. 
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Figure 12: Comparison of SUSY DGLAP and SUSY MC fragmentation functions for Mx = 
1 X 10^^ GeV with quark as a primary parton. SUSY DGLAP FFs are calculated for 10 TeV 
as the starting scale (solid line) and for Mz (broken line). SUSY MC FF is shown by dotted 
line. 
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Figure 13: Ratio of FFs calculated by MC and DGLAP in ordinary QCD (solid curve) and 
SUSY QCD (dashed curve) for Mx = 1 x 10^^ GeV. For discussion see the text. 
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Figure 14: Fraction of pions, e^, and nucleons, relative to all hadrons according to SUSY 
MC simulations. For pions the curve a corresponds to Mx = Mz, while curve b describes with 
a good accuracy the scales 10^° < Mx < 10^^ GeV. For nucleons curves c, d, e, / correspond to 
Mx equal to 10^^ 10^°, 10^ GeV and Mz, respectively 
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Figure 15: Neutrino (upper curves), gamma (middle curves) and nucleon (lower curves) spectra 
from DGLAP evolution at scales = 1 x lO^^ GeV (solid line) and Mx = 1 x 10^^ GeV 
(dashed line). The spectra are obtained averaging over primary partons. 




Figure 16: Comparison of photon spectra from present work computed with DGLAP equa- 
tion (solid line), from |33] (dashed line) and Pl] (dotted line). All three calculations are 
performed with quark as initial parton for Mx = 1 x 10^^ GeV. 
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Figure 17: Comparison of proton spectra from present work computed with MC (solid line), 
from 11211111 (dashed line) and from [Sll (dotted line). All spectra are computed with quark 
as initial parton for Mx = 1 x 10^^ GeV. 
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Figure 18: Gamma/nucleon ratio in generation spectra for Mx = 1 x 10^'* GeV computed with 
MC. The dashed curves illustrate the uncertainties of calculations. Calculations with DGLAP 
give very similar results. 
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Figure 19: Spectra of neutrinos (upper curve), photons (middle curve) and protons (two lower 
curves) in SHDM model compared with AGASA data. The neutrino flux is dominated by the 
halo component with small admixture of extragalactic flux. The flux of extragalactic protons 
is shown by the lower curve. 
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Figure 20: Comparison of SHDM prediction with the AG ASA data. The calculated spectrum 
of SHDM photons is shown by dotted curves for two different normalizations. The dashed 
curve gives the spectrum of extragalactic protons in the non-evolutionary model of Ref. jTO] . 
The sum of these two spectra is shown by the thick curves. The values are given of the 
comparison of these curves with experimental data for E' > 4 x 10^^ eV. 
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Figure 21: Diffuse spectra from necklaces. The upper curve shows neutrino flux,the middle 
- proton flux, and two lower curves - photon fluxes for two cases of absorption. The thick 
continuous curve gives the sum of the proton and higher photon flux. 



